/*
 * case.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#include "case_PAM_s1_h0_mu.h"
// MY HEADERS
#include"mros-v42.h"

namespace case_PAM_s1_h0_mu
{
global glob;

// *******************************************************************************************  //
// *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY ***
// *******************************************************************************************  //

#define Sqrt sqrt
#define Power pow

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double mu;
	double t;
	double tPrim;
	double U;
	double V;
	//double h;
	double epsilonF;
	double beta;

	unsigned int L;
};

struct variables {
	double lnF;
	//double ln;
	//double lmC;
	//double lmF;
	double nF;			// ok
	//double n;			// ok
	//double mC;			//
	//double mF;			//
	double dF;			// ok

    double Energy;
    double n, m;
 };

/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

inline void clear(parameters &p) {
	p.mu = p.t = p.tPrim = p.U = p.V = p.epsilonF = p.beta = 0.; p.L = 0;
}

inline void clear(variables &v) {
	v.lnF = v.nF = v.dF = 0;
	v.Energy = v.n = v.m = 0;
}

inline void addTo(const variables &q, variables &v, double d=1.) {
	v.lnF += q.lnF*d;
	v.nF += q.nF*d;
	v.dF += q.dF*d;
	v.Energy += q.Energy*d;
	v.n += q.n*d;
	v.m += q.m*d;
}

inline void copyTo(variables &v, const variables &q, double d=1.) {
	v.lnF = q.lnF*d;
	v.nF = q.nF*d;
	v.dF = q.dF*d;
	v.Energy = q.Energy*d;
	v.n = q.n*d;
	v.m = q.m*d;
}

inline void copyTo(parameters &v, const parameters &q) {
	v.mu = q.mu;
	v.t = q.t;
	v.tPrim = q.tPrim;
	v.U = q.U;
	v.V = q.V;
	v.epsilonF = q.epsilonF;
	v.beta = q.beta;
	v.L = q.L;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters  p.n  p.t  p.tPrim  p.U  p.V  p.h  p.epsilonF  p.beta  p.L
	variables v; // variables  v.lnF  v.lnF  v.lmC  v.lmF  v.nF  v.mu  v.mC  v.mF  v.d  v.Energy
	double nF, dF, lnF; // contractions

	// Abbreviations - dependend on kx, ky
	double epsilonk;

public:
	double E[4];
	double dEdlnF[4];
	double dEdq[4];

	double q[4];
	double dqdnF[4];
	double dqddF[4];

	Func(const parameters &p, const variables &v) {
		set(p, v);
		//refresh();
	}

	void set(const parameters &p, const variables &v) {
		copyTo(this->p, p);
		copyTo(this->v, v);
		refresh();
	}

	void set(const parameters &p) {
		copyTo(this->p, p);
		refresh();
	}

	// Gutzwiller
	double fq(int s);
	double fdqdnF(int s);
	double fdqddF(int s);

    // Energy
    double fE(double epsilonk, int s, int l);
    double fdEdlnF(double epsilonk, int s, int l);
    double fdEdq(double epsilonk, int s, int l);

	// REFRESH
	void refresh() {
		nF = v.nF;
		dF = v.dF;
		lnF = v.lnF;

		// obliczam funkcje Gutzwillara i poszczegolne pochodne
		for(int i=0; i<2; i++) {
			for(int j=0; j<2; j++) {
				q[j+2*i] 	 = fq(1-2*i);
				dqdnF[j+2*i] = fdqdnF(1-2*i);
				dqddF[j+2*i] = fdqddF(1-2*i);
			}
		}

		// zeby pozniej mozna bylo wypisac
		glob.qup =  q[0];
		glob.qupdnF =  dqdnF[0];
		glob.qupddF =  dqddF[0];

	}

	//REFRESH XY
	void refresh(double kx, double ky) {
		epsilonk = 2.*p.t*(cos(kx) + cos(ky)) + 4.*p.tPrim*cos(kx)*cos(ky);

		refresh();

		for(int i=0; i<2; i++) {
			for(int j=0; j<2; j++) {
				E[j+2*i] 	  = fE(epsilonk, 1-2*i, 1-2*j);
				dEdlnF[j+2*i] = fdEdlnF(epsilonk, 1-2*i, 1-2*j);
				dEdq[j+2*i]   = fdEdq(epsilonk, 1-2*i, 1-2*j);
			}
		}
	}
};


/////////////////
///////////////// FullSimplify
/////////////////

double Func::fq(int s) {
	double mF = 0;
	return (2.*Power(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		       Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))),2))/((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double Func::fdqdnF(int s) {
	double mF = 0;
	return ((dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		       Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
		     (4.*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((0.25*dF)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          (0.25 + 0.75*Power(dF,2) - 0.5*nF - 0.25*mF*s)/
		           Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) +
		       1.*(nF + mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
		       2.*(1 - 0.5*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
		   (Power(nF + mF*s,2)*Power(1 - 0.5*(nF + mF*s),2));
}

//double Func::fdqdmF(int s) {
//	return ((dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
//		       Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
//		     (4.*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((-0.25*dF*s)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
//		          (0.25*(1 + Power(dF,2) - nF)*s)/Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) +
//		       1.*s*(nF + mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
//		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
//		       2.*s*(1 - 0.5*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
//		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
//		   (Power(nF + mF*s,2)*Power(1 - 0.5*(nF + mF*s),2));
//}

double Func::fdqddF(int s) {
	double mF = 0;
	return 0. + (4.*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
	        Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
	      (0. - Power(dF,2)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
	        Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
	        (0. - 4*Power(dF,3) + dF*(-2. + 3.*nF + 1.*mF*s))/
	         (2.*Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
	    ((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double Func::fE(double epsilonk, int s, int l) {
	double ln = 0; double lmC = 0; double lmF = 0;
	return (epsilonk - 2*ln - lnF - 2*p.mu + p.epsilonF -
		     l*Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)))/2.;
}

double Func::fdEdlnF(double epsilonk, int s, int l) {
	double lmC = 0; double lmF = 0;
	return (-(l*(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s)) -
		     Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)))/
		   (2.*Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)));
}

double Func::fdEdq(double epsilonk, int s, int l) {
	double lmC = 0; double lmF = 0;
	return -((l*Power(p.V,2))/Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)));
}

/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	void onePoint(variables &y, int i, int j);

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->set(p);
	}

	void setParVar(parameters &p, variables &v) {
		this->p = p;
		this->v = v;
		f->set(p, v);
	}

	void run(variables &w);
};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	f->refresh(x, y);
	//if(DEEP_DEBUG) { f->printLambda(); }
	// DEBUG!!!
	if(i == 45 && j == 35) { glob.dEdq =  f->dEdq[0]; glob.dEdlnF =  f->dEdlnF[0];}

	double Etan[4] = { (1.-tanh(p.beta*f->E[0]/2.))/2., (1.-tanh(p.beta*f->E[1]/2.))/2., (1.-tanh(p.beta*f->E[2]/2.))/2., (1.-tanh(p.beta*f->E[3]/2.))/2.};

	for(int i=0; i<4; i++) {
		w.lnF += Etan[i]*f->dEdq[i]*f->dqdnF[i];
		w.nF += Etan[i]*f->dEdlnF[i];
		w.n += Etan[i];
		w.dF += Etan[i]*f->dEdq[i]*f->dqddF[i];

		if(i==0 || i==1) { w.m += Etan[i]; }
		else { w.m -= Etan[i]; }

        //if(f->E[i]>0) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
        //if(f->E[i]<0) { w.Energy += f->E[i]; }

		if(DEEP_DEBUG) { printf("########### %f, %f; %d: %f | %e, %e\n", x, y, i, Etan[i], f->E[i], p.beta*f->E[i]); }

        // it would be safer have following code - warning: 30 is chosen arbitrarily
		if(p.beta*f->E[i] > 30.) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
		else if(p.beta*f->E[i] < -30) { w.Energy += f->E[i]; }
		else { w.Energy += -log(1 + exp(-p.beta*f->E[i]))/p.beta; }
	}
	if(DEEP_DEBUG) {
		printf("###--- lnF=%.15e nF=%.15e dF=%.15e || %.15e %.15e %.15e %.15e | %.15e\n",
				w.lnF, w.nF, w.dF, f->E[0], p.beta*f->E[0], -exp(-p.beta*f->E[0])/p.beta, -log(1 + exp(-p.beta*f->E[0]))/p.beta, Etan[0]);
	}
}

void Equa::run(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation
	unsigned int max = p.L/2; // attention: L should be even!
	//glob.n = 0;

    // SMART WAY - ALSO OK
    // SERIAL CASE
   	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
   	// * edges (4 times)
   	for(unsigned int i=1; i<max; i++) {
   		onePoint(temp, i, i);
   		onePoint(temp, max, i);
   		onePoint(temp, i, 0);
   		//iii += 12;
   	}
      addTo(temp, w, 4);
      clear(temp);
      // * interior area (8 times) - SERIAL CASE
      for(unsigned int i=1; i<max; i++) {
   		for(unsigned int j=1; j<i; j++) {
   			onePoint(temp, i, j); //iii+=8;									//	sleep(1);
        	}
   	}
   	addTo(temp, w, 8);
   	clear(temp);
   	// * missing corners
   	onePoint(w, max, max);
   	onePoint(w, max, 0);
   	onePoint(w, 0, max);
   	onePoint(w, 0, 0);

	// what was out of sumation
	double L2 = (double)(p.L*p.L);
	w.lnF /= L2;
	w.nF /= L2;
	w.n /= L2;
	w.dF /= L2;

	w.Energy /= L2;
    w.m /= L2;

	w.lnF += v.lnF;
	w.nF += v.nF;
	w.dF += 2.*p.U*v.dF;

    v.Energy = w.Energy + v.lnF*v.nF + p.U*v.dF*v.dF;

	//if(w.chiAB != w.chiAB) { w.chiAB = 1000; }
	//if(w.chiS != w.chiS) { w.chiS = 1000; }
	//if(w.chiT != w.chiT) { w.chiT = 1000; }
	//if(w.barDelta != w.barDelta) { w.barDelta = 1000; }
	//if(w.simDelta != w.simDelta) { w.simDelta = 1000; }
	//if(w.mu != w.mu) { w.mu = 1000; }
	//if(w.mAF != w.mAF) { w.mAF = 1000; }

	if(DEEP_DEBUG) {
		printf("###+++ lnF=%.15e nF=%.15e dF=%.15e\n",
				w.lnF, w.nF, w.dF);
	}

    // results
	case_PAM_s1_h0_mu::glob.n = w.n;
	case_PAM_s1_h0_mu::glob.m = w.m;
    case_PAM_s1_h0_mu::glob.Energy = v.Energy + p.mu*glob.n;
	//glob.qt = f->fgt();

    //sleep(1);
}

// ***********************************************************************************  //
// *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY ***
// ***********************************************************************************  //


int function(const gsl_vector * x, void *params, gsl_vector * f) {

	struct parameters p;
	p.mu = ((Par *) params)->mu();
	p.t = ((Par *) params)->t();
	p.tPrim = ((Par *) params)->tPrim();
	p.U = ((Par *) params)->U();
	p.V = ((Par *) params)->V();
	p.epsilonF = ((Par *) params)->epsilonF();
	p.beta = ((Par *) params)->beta();
	p.L = ((Par *) params)->L();

	variables v;
	v.lnF = gsl_vector_get (x, 0);
	v.nF = gsl_vector_get (x, 1);
	v.dF = gsl_vector_get (x, 2);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2));
	}

    // It's start of the calculations

    Equa E(p, v);
    variables w;
    E.run(w);

    // It's end of the calculations

	gsl_vector_set (f, 0, w.lnF);
	gsl_vector_set (f, 1, w.nF);
	gsl_vector_set (f, 2, w.dF);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2));
		printf("###  %24.15e# %25.15e# %24.15e# | E = %24.15e n = %24.15e\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2), case_PAM_s1_h0_mu::glob.Energy, case_PAM_s1_h0_mu::glob.n);
		printf("\n");
	}

    return GSL_SUCCESS;
}

void ismaximum(Par &inpar, Var &invar) {
	printf("#"); inpar.print_verboten();
	Var initV = invar;

	// Solving
	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var exact(*res.get_V());
	printf("#Wynik: \n#");
	res.print_label();
	root.print();
	printf("\n");
	initV = exact;

	// EXTRA
	//Var exact = invar;
	//initV = exact;

	// Checking
	printf("#is it maximum?\n");
	printf("#          lnF                  nF                       d\n");

	struct parameters p;
	p.mu = inpar.mu();
	p.t = inpar.t();
	p.tPrim = inpar.tPrim();
	p.U = inpar.U();
	p.V = inpar.V();
	p.epsilonF = inpar.epsilonF();
	p.beta = inpar.beta();
	p.L = inpar.L();

	variables v; clear(v);
	Equa E(p, v);

	// EXTRA
	//v.lnF = initV.lnF();
	//v.lmF = initV.lmF();
	//v.nF = initV.nF();
	//v.mu = initV.mu();
	//v.mF = initV.mF();
	//v.dF = initV.dF();
	//variables w; //wynik
	//E.setParVar(p,v);
	//E.run2(w);

	double OptE = glob.Energy;
	double OptEnomu = glob.Energy - inpar.mu()*glob.n;
	while(1) {
		initV = exact; initV.shake(1e-4);
		v.lnF = initV.lnF();
		v.nF = initV.nF();
		v.dF = initV.dF();

		variables w; //wynik
		E.setParVar(p,v);
		E.run(w);
		if(abs(w.lnF) > 1e-6) { continue; }

		double globE = glob.Energy;
		double globEnomu = glob.Energy - glob.n*p.mu;
		double dif = globE-OptE;
		double difnomu = globEnomu-OptEnomu;

		initV.print();
		printf("%22.15f %22.15f %22.15f | %13e", glob.Energy, glob.n, glob.m, dif); if(dif<0) { printf(" Max!"); } else { printf("     "); }
		printf(" | (nomu): %13e || %e", difnomu, abs(w.lnF)); if(difnomu<0) { printf(" Max!"); }
		printf("\n");
		fflush(stdout);

		// EXTRA
		//if(difnomu < 0) {
		//	OptE = globE; OptEnomu = globEnomu;
		//	initV.set_lnF(v.lnF);
		//	initV.set_nF(v.nF);
		//	initV.set_dF(v.dF);
		//}
	}

	printf("#end ismaximum\n\n\n");
	fflush(stdout);
	return;
}

// ************* CLASSES ************* //

// Par

void Par::print(FILE *file) const {
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
}

void Par::print_verboten(FILE *file) const {
	fprintf(file, "PARAMETERS: mu=%f, t=%f, tPrim=%f, U=%f, V=%f, epsilonF=%f, beta=%f, L=%d\n", mu(), t(), tPrim(), U(), V(), epsilonF(), beta(), L());
}

void Par::print_label(FILE *file) const {
	fprintf(file, "mu   t   tPrim   U   V   epsilonF  beta   L\n");
}

// Var

void Var::print(FILE *file) const{
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
}

void Var::print_verboten(FILE *file) const{
	fprintf(file, "VARIABLES:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
	fprintf(file, "\n");
}

void Var::print_label(FILE *file) const {
	fprintf(file, "        lnF                  nF                     d");
}

// Res

void Res::set(gsl_multiroot_fsolver *s) {
	for(int i=0; i<n; i++) {
		(*V)[i] = gsl_vector_get (s->x, i);
	}
	this->glob = case_PAM_s1_h0_mu::glob;
}

void Res::print(FILE *file) const { V->print(file); fprintf(file, "%21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f\n", glob.Energy, glob.n, glob.m, glob.qup, glob.qupdnF, glob.qupddF, glob.dEdq, glob.dEdlnF); }
void Res::print_verboten(FILE *file) const { fprintf(file, "Result: "); V->print_verboten(file); fprintf(file, "%21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f\n", glob.Energy, glob.n, glob.m, glob.qup, glob.qupdnF, glob.qupddF, glob.dEdq, glob.dEdlnF); }
void Res::print_label(FILE *file) const { V->print_label(file); fprintf(file, "                      E                      n                      m                      qup                    qupdnF             qupddF                dEdq               dEdlnF\n");}

void example_inf() {
    Par par(N_PAR);
    par.set_mu(-1.769607670662955e+00);
    par.set_t(-1.);
    par.set_tPrim(0.);
    par.set_U(9);
    par.set_V(-2);
    par.set_epsilonF(-4);
    par.set_beta(1500.);
    par.set_L(256);
    par.print_verboten();

	// I'm choosing start point
	Var initV(N_VAR);

	FILE * pFile;
	pFile = fopen(OUTPUT_FILE, "w");

	// INFINITE
	Res res(N_VAR);
	RootGSL root(par, function, res);
	printf("sdelta      "); res.print_label(); fflush(stdout);
	while(1) {
	    initV.set_lnF(8*frand()-4);
	    initV.set_nF(2.*frand());
	    initV.set_dF(frand());

		if(root.Solve(initV) == GSL_SUCCESS) {
			//printf("       "); root.print(); root.print(pFile); par.print_verboten(); fflush(pFile); fflush(stdout);
			printf("       "); root.print(); root.print(pFile); fflush(pFile); fflush(stdout);

		}
	}
	fclose (pFile);
}

void run_1D() {

    // STANDARD
    Par initPar_st(N_PAR);
    initPar_st.set_mu(-1.769607670662955e+00);
    initPar_st.set_t(-1.);
    initPar_st.set_tPrim(0.);
    initPar_st.set_U(9);
    initPar_st.set_V(-2);
    initPar_st.set_epsilonF(-4.);
    initPar_st.set_beta(1500.);
    initPar_st.set_L(256);

//#automatic_run_n, start 1.600000, stop 2.000000, step 0.010000 | PARAMETERS: n=1.600000, t=-1.000000, tPrim=0.000000, U=9.000000, V=-2.000000, h=0.000000, epsilonF=-4.000000, beta=1500.000000, L=256
//#n                   lnF                  lmF                  nF                     mu                    mF                   d                      E                      m                      qup                    qdown
//    -3.118177359602460e+00  8.498622921750297e-01 8.807312115354573e-01 -1.489699194159749e+00  4.651329260873914e-01 2.224798529111916e-01   -6.466131210192415     0.400000000000000     0.774554817843601     0.698580678229792
//    -3.118177359602460e+00 -8.498622921750324e-01 8.807312115354576e-01 -1.489699194159749e+00 -4.651329260873924e-01 2.224798529111916e-01   -6.466131210192417    -0.400000000000001     0.698580678229792     0.774554817843601
//    -2.974984778852906e+00 -6.059331342162803e-17 8.673302395656413e-01 -1.769607670662955e+00 -4.782985343746993e-17 2.296590836132194e-01   -6.400003996300072     0.000000000000000     0.676132441880911     0.676132441880911

    Var initV_st(N_VAR);
    initV_st.set_lnF(-2.974984778852906e+00);
    initV_st.set_nF(8.673302395656413e-01);
    initV_st.set_dF(2.296590836132194e-01);

    Var initV_test(N_VAR);
    initV_test.set_lnF(-2.965888748423958e+00);
    initV_test.set_nF(8.577919263360242e-01);
    initV_test.set_dF(2.365572081178440e-01);

	// I'm choosing start point
    Par initPar(N_PAR); initPar = initPar_st;
	Var initV(N_VAR); initV = initV_test;

	ismaximum(initPar, initV);

	printf("#dochodze do punktu startowego\n");

	// STANDARD
	//automatic_run_n(initPar, initV, 1.75, 1.90, 0.01); printf("\n");
	//automatic_run_epsilonF(initPar, initV, -2, -5.0, 0.1); printf("\n");
	//automatic_run_h(initPar, initV, 0.03, 0.002, 0.002); printf("\n");
	//automatic_run_V(initPar, initV, -3, -2.0, 0.2); printf("\n");
	//automatic_run_tPrim(initPar, initV, 0.25, 0.0, 0.05); printf("\n");
	//automatic_run_U(initPar, initV, 0, 4., 1); printf("\n");
	//automatic_run_U(initPar, initV, 1., 0., 0.1); printf("\n");
	//automatic_run_U(initPar, initV, 0.1, 0.01, 0.01); printf("\n");
	//automatic_run_n(initPar, initV, 1.75, 2.50, 0.01); printf("\n");
	//automatic_run_U(initPar, initV, 1, 12.0, 0.01); printf("\n");

	printf("#Jestem w punkcie startowym\n");
    Par oldPar = initPar;
    Var oldVar = initV;

    // STANDARD
    //automatic_run_n(initPar, initV, 2.00, 0.00, 0.05); initPar=oldPar; initV=oldVar;


    // ZERO
    automatic_run_mu(initPar, initV, -1.769607670662955e+00, -6, 0.005); initPar=oldPar; initV=oldVar;


	return;
}

void automatic_run_mu(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_n, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#mu           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_mu(start + ((double)i)*step);
		printf("%f ", inpar.mu());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	inpar.set_mu(stop);
	invar = initV;

	printf("#end automatic_run_n\n\n\n");
	fflush(stdout);
	return;
}

void automatic_run_tPrim(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_tPrim, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#tPrim       "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_tPrim(start + ((double)i)*step);
		printf("%f ", inpar.tPrim());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_tPrim\n\n\n");
	fflush(stdout);

	inpar.set_tPrim(stop);
	invar = initV;
	return;
}

void automatic_run_U(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_U, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#U           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_U(start + ((double)i)*step);
		printf("%f ", inpar.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_U\n\n\n");
	fflush(stdout);

	inpar.set_U(stop);
	invar = initV;
	return;
}

void automatic_run_V(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_V, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#V           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_V(start + ((double)i)*step);
		printf("%f ", inpar.V());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_V\n\n\n");
	fflush(stdout);

	inpar.set_V(stop);
	invar = initV;
	return;
}

void automatic_run_epsilonF(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_epsilonF, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#epsilonF    "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_epsilonF(start + ((double)i)*step);
		printf("%f ", inpar.epsilonF());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_epsilonF\n\n\n");
	fflush(stdout);

	inpar.set_epsilonF(stop);
	invar = initV;
	return;
}

void automatic_run_beta(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_beta, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#beta        "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_beta(start + ((double)i)*step);
		printf("%f ", inpar.beta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_beta\n\n\n");
	fflush(stdout);

	inpar.set_beta(stop);
	invar = initV;
	return;
}

void automatic_run_L(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_L, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#L           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_L(start + ((double)i)*step);
		printf("%d ", inpar.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_L\n\n\n");
	fflush(stdout);

	inpar.set_L(stop);
	invar = initV;
	return;
}

} // end of namespace case2D2t_SGA_tPrim0
